y5 = y5 - y5[ind] + vertical; y6 = y6 - y6[ind] + vertical; y7 = y7 - y7[ind] + vertical;
y8 = y8 - y8[ind] + vertical;
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_birthweight_OLS_poly.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2900,3300),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2900, 50, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2900,3300,50),labels =seq(2900,3300,50), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y2, col = 'red3', lwd=2, lty=1) #quadratic
lines(x, y3, col = 'gold', lwd=2, lty=1) #cubic
lines(x, y4, col = 'forestgreen', lwd=2, lty=1) #quartic
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_birthweight_OLS_FE.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2900,3300),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2900, 50, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2900,3300,50),labels =seq(2900,3300,50), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y5, col = 'red3', lwd=2, lty=1) #tpu_year no controls FE
lines(x, y6, col = 'gold', lwd=2, lty=1) #tpu_year controls FE
lines(x, y7, col = 'forestgreen', lwd=2, lty=1) #tpu_year month FE
lines(x, y8, col = 'purple', lwd=2, lty=1) #mother FE
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
rm(list=ls())
source("/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/AnalysisCode/helper_functions.R")
data <- read.csv("/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/ProcessedData/Final/Mortality_NonUniqueMatch_10km_samples_stored.csv")
#Specify model components
regression_controls <- "maxtemp + mintemp + rain + meanhum + pressure"
regression_fixed_effects <- "tpu_year + tpu_month + MomAge + MomEdu + MomActive + MomMaritalStatus + DadActive + DadAge + HousingType + Girl + TypeBirth + NumPreviousBirth + Hospital2"
regression_var_cluster <- "TPU"
summary(data$RSP)
quantile(data$RSP, c(.90, .99), na.rm=TRUE)
pm10_support <- 30:80
main_model_results_a <- lfe::felm(BirthWeight ~ maxtemp + mintemp + rain + meanhum + pressure | tpu_year + tpu_month + MomAge + MomEdu + MomActive + MomMaritalStatus + DadActive + DadAge + HousingType + Girl + TypeBirth + NumPreviousBirth + Hospital2 | (RSP ~ major_inversion) | TPU + date_of_birth, data = data, Nboot = 1000, bootexpr = quote(summary(est)$coefficients["`RSP(fit)`","Estimate"]*(pm10_support)))
bstrap_a <- main_model_results_a$boot #pull out bootstrap runs
#use bootstrapped model runs to generate 95% confidence intervals for imr | pm  across inner 99 percentile of observed pm range
bstrap_a <- apply(bstrap_a, 2, function(x){center(x,pm10_support, round(mean(data$RSP)), mean(data$BirthWeight))}) #re-center to (mean pm, mean imr)
response_a <- tibble::tibble(x =pm10_support,  y = apply(bstrap_a, 1, median), #bring in x and y
ci_low = apply(bstrap_a, 1, function(x){quantile(x, 0.025)}), ci_high = apply(bstrap_a, 1, function(x){quantile(x, 0.975)}) ) #%>%  #generate ci
#filter(x > quantile(data$RSP, 0.01,na.rm = TRUE) & x < quantile(data$RSP, 0.99, na.rm=TRUE))  #keep inner 99 percentile observed pm concentrations
filter(x >= 30 & x <= 80)  #keep inner 99 percentile observed pm concentrations
#re-organize x and y value into structure for plotting with polygon()
px <- c(response_a$x[1], response_a$x, response_a$x[nrow(response_a):1])
py <- c(response_a$ci_low[1], response_a$ci_high, response_a$ci_low[nrow(response_a):1])
ci_polygon_a <- tibble::tibble(px, py)
#define object that stores info on breaks and heights of histogram at bottom of panel (a)
hist_ht_a <- get_hst_obj(data$RSP,30,80,0.5)
######## define plot colors #########
#Derived from Paul Tol's palettes designed to be distinct in colour-blind vision
#https://personal.sron.nl/~pault/
col_line_a <- "#114477"
col_shade_a <- "#77AADD"
col_line_b1 <- "#114477"
col_shade_b1 <- "#77AADD"
col_line_b2 <- "#117744"
col_shade_b2 <- "#44AA77"
col_line_c1 <-"#DDAA77"
col_shade_c1 <-"#AA7744"
col_line_c2 <- "#AA4488"
col_shade_c2 <- "#CC99BB"
col_d <- "black"
col_e <- "black"
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_BirthWeight_IV.pdf",width = 9, height= 9, useDingbats = F)
x <- 30:80
###############################
#### Panel a: Main Effect
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2900,3300),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 3, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2900, 100, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2900,3300,100),labels =seq(2900,3300,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
abline(h=mean(data$BirthWeight), xlim=c(30,80),lty=2)
dev.off()
save(ci_polygon_a, response_a, hist_ht_a, file="/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/ProcessedData/Final/Figure_response_RSP_BirthWeight_IV.RData")
load("/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/ProcessedData/Final/Figure_response_RSP_BirthWeight_IV.RData")
first_stage_8 <- lfe::felm(RSP ~ lnNRT + Dist_lnNRT + maxtemp + mintemp + rain + meanhum + pressure | tpu_year + tpu_month + MomAge + MomEdu + MomActive + MomMaritalStatus + DadActive + DadAge + HousingType + Girl + TypeBirth + NumPreviousBirth + Hospital2 |0| TPU + date_of_birth, data = data)
RSP_fit8 <- first_stage_8$fitted.values
second_stage_8 <- summary(felm(as.formula(paste("BirthWeight ~ RSP_fit8 +",
regression_controls, "|", "tpu_year + tpu_month + MomAge + MomEdu + MomActive + MomMaritalStatus + DadActive + DadAge + HousingType + Girl + TypeBirth + NumPreviousBirth + Hospital2", sep = "")), data = data))
x <- response_a$x
ind <-which(x==round(mean(data$RSP,na.rm=TRUE)))
y8 = c((second_stage_8)$coefficients[1:1,"Estimate"])*x
vertical <-mean(data$BirthWeight)
y8 = y8 - y8[ind] + vertical;
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_BirthWeight_IV_both_IVs.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2900,3300),xlim = c(30,),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2900, 100, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2900,3300,100),labels =seq(2900,3300,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y8, col = 'red3', lwd=2, lty=1) #port traffic IV
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_BirthWeight_IV_both_IVs.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2900,3300),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2900, 5-, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2900,3300,100),labels =seq(2900,3300,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y8, col = 'red3', lwd=2, lty=1) #port traffic IV
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2600,3600),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2600, 100, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2600,3600,100),labels =seq(2600,3600,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y8, col = 'red3', lwd=2, lty=1) #port traffic IV
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_BirthWeight_IV_both_IVs.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2600,3600),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2600, 100, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2600,3600,100),labels =seq(2600,3600,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y8, col = 'red3', lwd=2, lty=1) #port traffic IV
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
pdf(file = "/Users/jonathancolmer/Dropbox/AirQuality&BabyBirth/Submission/Figures/PM10_BirthWeight_IV_both_IVs.pdf",width = 9, height= 9, useDingbats = F)
###############################
#### Different Specifications
################################
par(mgp=c(2,0.8,0),mar = c(5,4,4,3), lend=1)
#plot
plot(response_a$x, response_a$y,type = 'l',col = NA, lwd = 3,ylim = c(2700,3500),xlim = c(30,80),
axes = F,xlab = "",ylab = "")
polygon(ci_polygon_a$px, ci_polygon_a$py, col = add.alpha(col_shade_a,0.25),border =NA)
lines(response_a$x, response_a$y, lwd = 2, col = col_line_a)
plotHist(hist_ht_a,col_shade_a, 0.25, 2700, 100, add.alpha(col_line_a, 0.15) )
#axis
axis(1, tick=T, at = seq(30,80,5),col.axis='gray30',line=-1,cex.axis = 1.25)
axis(2, tick = T, at = seq(2700,3500,100),labels =seq(2700,3500,100), las =2 ,line = -1, cex.axis=1.25)
#labels
mtext(side = 1, text = "PM10 Concentration (ug/m3)",line = 1.5,cex = 1.5)
mtext(side = 2, text = "Birthweight (grams)",line = 2.5,cex = 1.5)
lines(x, y8, col = 'red3', lwd=2, lty=1) #port traffic IV
abline(h=mean(data$BirthWeight), lty=2)
dev.off()
username <- "/Users/jonathancolmer/"
directory <- "The Lab Dropbox/Jonathan Colmer/JMP_Submission/Data/ASI/"
read.csv("bunching_data.csv")
username <- "/Users/jonathancolmer/"
directory <- "The Lab Dropbox/Jonathan Colmer/JMP_Submission/Data/ASI/Analysis/"
read.csv("bunching_data.csv")
read.csv("/Users/jonathancolmer/The Lab Dropbox/Jonathan Colmer/JMP_Submission/Data/ASI/Analysis/bunching_data.csv")
data <- read.csv("/Users/jonathancolmer/The Lab Dropbox/Jonathan Colmer/JMP_Submission/Data/ASI/Analysis/bunching_data.csv")
install.packages("bunchr")
library(bunchr)
bunchApp()
install.packages('devtools')
library(devtools)
install_github('andreacirilloac/updateR')
library(updateR)
updater(admin_passworkd=Remloc89)
updateR(admin_passworkd=Remloc89)
updateR(admin_password=Remloc89)
updateR(admin_password='Remloc89')
install.packages(as.vector(needed_packages))
debugSource('~/The Lab Dropbox/Jonathan Colmer/JMP_Submission/AEJ Submission/Code/Processing/Weather/UDEL/Calculate Zonal Statistics.R')
devtools::install_github('isciences/exactextractr')
library(sp)
library(raster)
library(sf)
library(exactextractr)
library(rgdal)
library(gtools)
library(foreign)
library(gdalUtils)
library(ncdf4)
username <- "/Users/jonathancolmer/"
directory <- "The Lab Dropbox/Jonathan Colmer/JMP_Submission/AEJ Submission/Data/"
path <- paste0(username, "/", directory,"/")
### GROUND STATION TEMPERATURE ### District - Month - Year ###
### Read in list of ground monitor temperature data files ###
temp_years <- list.files(paste0(path,"/Data - Raw/air_temp_2017"),"*.", full.names = FALSE)
### We start with 51 because it refers to 1950 in the list of years
###GS Temp Loop###
for(i in 51:length(temp_years)){
india <- shapefile(paste0(path,"Data - Raw/Shapefiles/District_Boundaries.shp"))
# Keep only necessary attributes
india <- india[c(1, 3, 15:17, 22, 28)]
head(india)
### Get District Centroids ###
india$coords <- coordinates(india)
india$lon <- india$coords[,1]
india$lat <- india$coords[,2]
india$coords <- NULL
### Create sf object from the India Shapefile ###
india_sf <- st_as_sf(india)
# Read in the first air temp dataframe
df <- read.table(paste0(path,"/Data - Raw/air_temp_2017/",temp_years[i]),header = FALSE)
### These files are organized such that:
### columns 1 & 2 contain lon and lat,
### 3 - 14 contain monthly measures,
### and 15 contains annual mean temp.
df <- df[1:14]
colnames(df) <- c("lon", "lat", "temp_jan", "temp_feb", "temp_mar", "temp_apr","temp_may", "temp_jun", "temp_jul", "temp_aug", "temp_sep", "temp_oct", "temp_nov", "temp_dec")
### Generate list of column names to rasterize ###
temp_list <- colnames(df)
### Keep only the month labels in the list
temp_list <- temp_list[3:14]
### Add year column and assign year from file name
india$temp_year <- as.numeric(paste0(gsub("^.*\\.","",temp_years[i])))
### Create a raster object ###
resolution <- 0.5
xmin <- min(df$lon) - (.5 * resolution)
xmax <- max(df$lon) + (.5 * resolution)
ymin <- min(df$lat) - (.5 * resolution)
ymax <- max(df$lat) + (.5 * resolution)
temp_raster <- raster(xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax, res = resolution, crs=crs(india))
### Rasterize each column of monthly temp values and run zonal statistics with exactextractr package ###
for(j in 1:12){
temp_data <- rasterize(df[, c('lon', 'lat')], temp_raster, df[, paste0(temp_list[j])], fun=mean)
india$temp <- exact_extract(temp_data, india_sf, 'mean')
colnames(india@data)[colnames(india@data) == "temp"] <- paste0(temp_list[j])
}
### Convert to dataframe
data_frame <- as.data.frame(india)
setwd(paste0(path,"/Data - Clean/GS Temp by Year"))
### write data to yearly stata files
write.dta(data_frame, paste0(gsub("^.*\\.","",temp_years[i]),".dta"))
}
### GROUND STATION PRECIPITATION ### District - Month - Year ###
rain_years <- list.files(paste0(path,"/Data - Raw/precip_2017"),"*.", full.names = FALSE)
###GS Rain Loop###
for(i in 49:length(rain_years)){
india <- shapefile(paste0(path,"Data - Raw/Shapefiles/District_Boundaries.shp"))
india <- india[c(1, 3, 15:17, 22, 28)]
head(india)
### Get District Centroids ###
india$coords <- coordinates(india)
india$lon <- india$coords[,1]
india$lat <- india$coords[,2]
india$coords <- NULL
head(india)
### Create sf object from the India Shapefile ###
india_sf <- st_as_sf(india)
df <- read.table(paste0(path,"/Data - Raw/precip_2017/",rain_years[i]),header = FALSE)
### These files are organized such that:
### columns 1 & 2 contain lon and lat,
### 3 - 14 contain monthly measures,
### and 15 contains annual mean precipitation.
df <- df[1:14]
colnames(df) <- c("lon", "lat", "rain_jan", "rain_feb", "rain_mar", "rain_apr","rain_may", "rain_jun", "rain_jul", "rain_aug", "rain_sep", "rain_oct", "rain_nov", "rain_dec")
rain_list <- colnames(df)
rain_list <- rain_list[3:14]
india$rain_year <- as.numeric(paste0(gsub("^.*\\.","",rain_years[i])))
### Create a raster object ###
resolution <- 0.5
xmin <- min(df$lon) - (.5 * resolution)
xmax <- max(df$lon) + (.5 * resolution)
ymin <- min(df$lat) - (.5 * resolution)
ymax <- max(df$lat) + (.5 * resolution)
rain_raster <- raster(xmn = xmin, xmx = xmax, ymn = ymin, ymx = ymax, res = resolution, crs=crs(india))
for(j in 1:12){
rain_data <- rasterize(df[, c('lon', 'lat')], rain_raster, df[, paste0(rain_list[j])], fun=mean)
india$rain <- exact_extract(rain_data, india_sf, 'mean')
colnames(india@data)[colnames(india@data) == "rain"] <- paste0(rain_list[j])
}
data_frame <- as.data.frame(india)
setwd(paste0(path,"/Data - Clean/GS Rain by Year"))
### write data to yearly stata files
write.dta(data_frame, paste0(gsub("^.*\\.","",rain_years[i]),".dta"))
}
source('~/The Lab Dropbox/Jonathan Colmer/JMP_Submission/AEJ Submission/Code/Processing/Weather/UDEL/Calculate Zonal Statistics.R')
devtools::install_github('isciences/exactextractr')
install.packages("geos-config", type="source")
install.packages('rgeos')
install.packages('rgdal')
devtools::install_github('isciences/exactextractr')
install.packages('geos-config')
load("/Users/jonathancolmer/Downloads/PM25_1997_2016.rda")
write.csv(pm25_allyear,"/Users/jonathancolmer/Downloads/PM25_1997_2016.csv")
source('~/.active-rstudio-document')
source('~/.active-rstudio-document')
source('~/.active-rstudio-document')
library("dplyr")
library("pscl")
library("MASS")
library(NBZIMM)
library("lme4")
install.packages("NBZIMM")
install.packages("pscl")
glmm.zinb.off = glmm.zinb(fixed = Deaths ~ mean_pm25 +scale(poverty) + scale(popdensity)  +scale(medianhousevalue)
+scale(medhouseholdincome) + scale(pct_owner_occ)  +scale(hispanic)
+scale(education)  +scale(pct_blk) + scale(older_pecent)
+ scale(totalTestResults) +
+ scale(beds)
+ scale(mean_bmi) + scale(smoke_rate)
+ scale(mean_summer_temp) + scale(mean_winter_temp) + scale(mean_summer_rm) + scale(mean_winter_rm)
+ offset(log(population)),
random = ~ 1 | state, data = (aggregate_pm_census_cdc_test_beds))
library("dplyr")
library("pscl")
library("MASS")
library(NBZIMM)
library("lme4")
library("NBZIMM")
install.packages("NBZIMM")
install.packages(NBZIMM)
View(aggregate_pm_census)
View(aggregate_pm_census_cdc_test)
View(aggregate_pm_census_cdc_test_beds)
View(aggregate_pm_census_cdc_test_beds)
write.dta(aggregate_pm_census_cdc_test_beds, "Downloads/⁩")
write.dta(aggregate_pm_census_cdc_test_beds, "/Users/jonathancolmer/Downloads")
library("foreign")
write.dta(aggregate_pm_census_cdc_test_beds, "/Users/jonathancolmer/Downloads")
write.dta(aggregate_pm_census_cdc_test_beds, "/Users/jonathancolmer/Downloads/")
write.dta(aggregate_pm_census_cdc_test_beds, "/Users/jonathancolmer/Downloads/aggregate_pm_census_cdc_test_beds.dta")
libPaths()
libPaths(raster)
.libPaths(raster)
.libPaths()
print(plot3)
library(ggplot2)
# Generate some data
set.seed(123)
data1 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 60, max = 140))  # values between 60 and 140
data2 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 20, max = 80))   # values between 20 and 80
data3 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 100, max = 130)) # values between 100 and 130
# Create a function to generate the plots
generate_plot <- function(data) {
ggplot(data, aes(x = x, y = y)) +
geom_tile(aes(fill = z)) +
scale_fill_distiller(palette = "viridis",  # use viridis color palette
limits = c(0, 140))   # set limits of color scale to 0-140
}
# Generate the plots
plot1 <- generate_plot(data1)
plot2 <- generate_plot(data2)
plot3 <- generate_plot(data3)
# Print the plots
print(plot1)
print(plot2)
print(plot3)
library(virids)
library(viridis)
set.seed(123)
data1 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 60, max = 140))  # values between 60 and 140
data2 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 20, max = 80))   # values between 20 and 80
data3 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 100, max = 130)) # values between 100 and 130
# Create a function to generate the plots
generate_plot <- function(data) {
ggplot(data, aes(x = x, y = y)) +
geom_tile(aes(fill = z)) +
scale_fill_distiller(palette = "viridis",  # use viridis color palette
limits = c(0, 140))   # set limits of color scale to 0-140
}
# Generate the plots
plot1 <- generate_plot(data1)
plot2 <- generate_plot(data2)
plot3 <- generate_plot(data3)
# Print the plots
print(plot1)
print(plot2)
print(plot3)
install.packages("viridis")
library(viridis)
install.packages("viridis")
library(viridis)
# Generate some data
set.seed(123)
data1 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 60, max = 140))  # values between 60 and 140
data2 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 20, max = 80))   # values between 20 and 80
data3 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 100, max = 130)) # values between 100 and 130
# Create a function to generate the plots
generate_plot <- function(data) {
ggplot(data, aes(x = x, y = y)) +
geom_tile(aes(fill = z)) +
scale_fill_distiller(palette = "viridis",  # use viridis color palette
limits = c(0, 140))   # set limits of color scale to 0-140
}
# Generate the plots
plot1 <- generate_plot(data1)
plot2 <- generate_plot(data2)
plot3 <- generate_plot(data3)
# Print the plots
print(plot1)
print(plot2)
print(plot3)
install.packages("ggplot2")
install.packages("ggplot2")
install.packages("ggplot2")
install.packages("ggplot2")
install.packages("ggplot2")
# Generate some data
set.seed(123)
data1 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 60, max = 140))  # values between 60 and 140
data2 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 20, max = 80))   # values between 20 and 80
data3 <- data.frame(x = rnorm(100), y = rnorm(100), z = runif(100, min = 100, max = 130)) # values between 100 and 130
# Create a function to generate the plots
generate_plot <- function(data) {
ggplot(data, aes(x = x, y = y)) +
geom_tile(aes(fill = z)) +
scale_fill_distiller(palette = "viridis",  # use viridis color palette
limits = c(0, 140))   # set limits of color scale to 0-140
}
# Generate the plots
plot1 <- generate_plot(data1)
plot2 <- generate_plot(data2)
plot3 <- generate_plot(data3)
# Print the plots
print(plot1)
print(plot2)
print(plot3)
generate_plot <- function(data) {
ggplot2(data, aes(x = x, y = y)) +
geom_tile(aes(fill = z)) +
scale_fill_distiller(palette = "viridis",  # use viridis color palette
limits = c(0, 140))   # set limits of color scale to 0-140
}
# Generate the plots
plot1 <- generate_plot(data1)
plot2 <- generate_plot(data2)
plot3 <- generate_plot(data3)
# Print the plots
print(plot1)
print(plot2)
print(plot3)
library(ggplot2)
### PACKAGE PREAMBLE ###
# Create a vector of package names
packages <- c("dplyr", "ggplot2", "tidyr", "tidyverse", "data.table", "readxl", "foreign", "devtools", "fixest", "sandwich", "modelsummary", "kableExtra")
installed_packages <- packages %in% rownames(installed.packages())
if (any(installed_packages == FALSE)) {
install.packages(packages[!installed_packages])
}
# Packages loading
invisible(lapply(packages, library, character.only = TRUE))
mydata_UCR <- read.dta("UCR_analysis.dta")
setwd(dir = "/Users/jmc4qg/The Lab Dropbox/Jonathan Colmer/ShotSpotter_env/Journal_submissions/ReStat/Replication Folder/Analysis Data/")
setwd(dir = "/Users/jonathancolmer/The Lab Dropbox/Jonathan Colmer/ShotSpotter_env/Journal_submissions/ReStat/Replication Folder/Analysis Data/")
mydata_UCR <- read.dta("UCR_analysis.dta")
DiDiT_2_UCR <- feols(homicide_pc~ temp_MP + prec_MP | ori_year+ORI[tMean]+year[tMean]+ORI[prec]+year[prec], data=mydata_UCR)
install.packages("fixest")
install.packages("feols")
